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Abstract 


A  rigorous  formulation  of  rough  earth  parabolic  equation  (PE)  is  given.  The  formulation 
IS  based,  on  casting  the  governing  transform  equations  in  terms  of  incident  and  reflected 
plane  waves  and  using  the  rough-surface  reduction  factor  directly  in  the  spectral  do¬ 
main.  Solution  is  performed  by  the  well  known  split-step  algorithm.  Inclusion  of  surface 
roughness  into  the  PE  in  this  manner  requires  a  redefinition  of  the  Fourier  transform 
pair.  Several  examples  are  considered  showing  comparisons  with  waveguide  and  other 
methods. 


Objective 


Our  research  objective  is  to  develop  a  rigorous  method  to  incorporate  sea-surface  rough¬ 
ness  into  the  parabolic  equation  modeling  of  radiowave  propagation. 


The  Parabolic  Equation  (PE)  method  has  emerged  as  an  extremely  valuable  method 
for  assessing  radiowave  propagation  in  the  lower  atmosphere  in  the  presence  of  ducts. 
Propagation  loss  can  be  easily  estimated  over  very  long  ranges  of  the  order  of  a  few 
hundred  kilometers  for  frequencies  through  Super  High  Frequency  (SHF)  band,  and  for 
antenna  heights  extending  up  to  a  few  hundred  meters.  It  is  also  possible  to  directly 
account  for  finite  conductivity  of  the  earth  in  the  PE.  Reference  [1]  discusses  the  basic 
idea  and  the  various  approximations  involved  in  the  development  of  the  PE.  Although 
there  are  many  different  ways  of  solving  the  PE  [2],  none  seem  to  offer  the  computational 
advantages  of  the  split-step  algorithm,  originally  developed  by  Tappert  [3],  in  terms  of 
the  large  range  steps  allowable.  One  of  the  current  unresolved  issues  is  the  incorporation 
of  sea-surface  roughness  into  the  PE.  Several  attempts  have  been  made  in  the  past  all 
of  which  rely  on  some  kind  of  approximation  or  post-processing  of  data  obtained  by 
running  the  PE  over  smooth  earth.  Reference  [4]  discusses  some  of  these  methods  in 
detail.  In  this  report,  we  present  a  direct  way  of  incorporating  surface  roughness  into 
the  PE.  The  idea  behind  our  approach  is  to  cast  the  governing  transform  equations  in 
terms  of  incident  and  reflected  plane  waves,  and  then  use  the  rough  surface  reduction 
factor  available  for  plane  waves  [5]  directly  in  the  spectral  domain.  Of  course  this  will 
nesseciate  the  modification  of  the  Fourier  transform  pair,  as  we  will  show  shortly. 
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Theory 

1.  Smooth  Earth 


The  starting  point  for  our  formulation  is  the  standard  parabolic  equation,  given,  for 
example,  in  [1].  Assuming  an  time  dependence,  we  consider  the  parabolic  equation 
in  a  medium  with  parameters 

—  (x,  z)  +  2iA:o  — (x,  z)  +  2klM{x,  z)u{x,  z)  =  0,  (1) 

where  u(x,z)  =  V r  sin  6E^[r,  9),  for  horizontal  polarization  and  u(x,z)  =  \/q*r^) 
for  vertical  polarization,  (r,  9,  (j))  being  the  usual  spherical  coordinates.  The  coordinate 
system  is  chosen  such  that  the  source  is  located  at  9  =  0  and  r  =  +  z,,  where  is 

the  radius  of  the  earth  and  x,  is  the  height  of  the  source  relative  to  the  surface  of  the 
earth.  Furthermore,  ko  =  is  the  free-space  wavenumber,  x  is  the  range  axis,  z  is 

the  height  axis,  and 

M{x,  z)  =  —  1  +  — ^  X  10® 

is  the  modified  refractive  index  under  earth-flattened  approximations  [1].  The  quantity 
(H^)  is  the  (^-component  of  the  electric  (magnetic)  field.  The  complex  dielectric 
constant,  e^c,  of  the  earth  is  ere  =  Cr  +  icr/wco,  where  a  is  its  conductivity.  The  parabolic 
equation  given  in  (1)  is  to  be  solved  subject  to  the  approximate  boundary  condition 

^(x,  0)  -t-  aou{x,  0)  =  0,  (2) 


where  ao  =  iko^/e^7-^  for  horizontal  polarization,  ao  =  ikoy/e^c  —  l/^rc  for  vertical 
polarization.  The  above  boundary  condition  is  approximate  because  it  is  valid  for  highly 
conductive  soil  satisfying  Icrd  >  1.  If  this  condition  is  not  met,  equation  (2)  must  be 
imposed  seperately  for  each  constituent  plane  wave  comprising  u  with  oto  replaced  by 


a  =  iko  cos  9i 


1  -  r(erc;  g.) 
1  -b  r(erc;  9i) 


where  F  is  the  reflection  coefficient  for  a  plane  wave  striking  the  earth  at  an  angle  9i  with 
the  normal.  The  actual. field,  u,  does  not  satisfy  the  simple  condition  dictated  by  (2). 
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To  solve  (1)  subject  the  boundary  condition  (2),  the  following  mixed  Fourier  transform 
pair  is  defined  [1]: 


uix 


oo 

^p)  =  J^{u)  =  j  u(x,  z)[Q:osinp2  —  pcos  pz]  dz 


/  X  2  /•.,  ^aosinpz -pcospz 

u(x,z)  =  z  («)  =  -y»(x,p) — 


■  dp  +  S{x)e' 


(за) 

(зб) 


where  the  latter  term  in  (3b)  decays  with  height  and  range  and  is  due  to  the  surface 
wave  propagating  on  a  non-perfectly  conducting  earth.  This  term  is  usually  ignored 
in  the  solution  for  frequencies  >  10  MHz.  In  the  visible  range  —ko  <  p  <  one  may 
identify  p  with  the  vertical  wavenumber  ko  sin  V’,  where  ip  is  the  grazing  angle  with  respect 
to  horizontal  taken  positive  for  waves  approaching  the  interface  2  =  0. 

Using  the  above  Fourier  transform  pair,  the  solution  to  (1)  in  the  presence  of  a  slowly 
varying  duct  can  be  written  as 

u{x,  z)  =  e‘*o(x-ro)io-«xAf^-i  |e-*P^(^-^o)/2fcojr[u(a;o,  2)]}  (4) 

To  appreciate  the  development  of  the  rough  earth  case,  we  will  first  rewrite  the  above 
Fourier  transform  pair  (ignoring  the  surface  wave)  in  the  following  equivalent  form: 


00  I  ^ 

u{x,p)  =  =  J  u{x,z)e'^^  dz  +  j  uix,z)e-'‘ 

ti(z,2)  =  =  —  J  u{x,p)  +  r,(p)e’'’^]  dp, 


e-'P^  dz 

(5a) 

dp, 

(56) 

where 


r,(p)  =  ^i^,  p>o 

p-tao 


is  the  plane- wave  reflection  coefficient  for  smooth  earth,  and  subscripts  s  denote  smooth 
earth  case.  In  this  form,  it  is  clear  from  looking  at  (5b)  that  the  total  fleld  is  comprised 
of  an  incident  plane  wave  of  spectral  amplitude  u{x,p)  traveling  towards  the  interface 
2  =  0  and  a  corresponding  reflected  plane  wave  of  amplitude  Ts{p)u[x,p)  traveling  away 
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from  the  interface.  Furthermore,  the  field  u  satisfies  the  required  smooth  earth  boundary 
condition  of  (2).  Although  the  expression  given  above  for  Tslp)  makes  physical  sense 
only  for  p  >  0,  we  may  analytically  continue  it  for  p  <  0  and  use  the  same  expression  to 
define  it  for  negative  p.  With  this  continuation,  we  see  that  rj(— p)  =  l/rs(p).  A  second 
equivalent  form  which  is  symmetric  with  respect  to  the  reflection  coefficient  and  which 
makes  use  of  the  extended  definition  of  rs(p)  can  be  written  down  as: 

OO  CO 

u(x,p)  =  =  yrTrt  J  u(x,zy-dz--^  J  u(x,z)e-‘''-dz  (6a) 

0  sip)  0 

u(x,  z)  =  ^  ‘‘p-  («'’) 

Note  that  there  are  no  approximations  involved,  and  all  the  three  forms  of  the  Fourier 
transform  pair  are  equivalent.  Also  note  that  the  form  in  (6)  is  particularly  useful  in 
extending  the  PE  methodology  to  rough  surfaces.  The  transform  u  is  seen  to  be  an  odd 
function  of  the  transform  variable  p. 

2.  Rough  Earth 

In  the  rough  earth  case,  we  simply  replace  the  smooth  earth  spectral  reflection  coefficient 
V,{p)  with  the  corresponding  rough  earth  one  in  the  inversion  formula  (6b).  A  forward 
transform  formula  consistent  with  this  must  be  defined  seperately.  The  rough  surface 
reflection  coefficient,  rr(p),  is  taken  as 

^r{p)  =  Po{p]<^h)Ts{p),  (^) 

where  Po(p;<7’^)  is  the  rough  surface  reduction  factor  defined  by  [5] 

Po(p;^/.)  =  e-2pM/o(2pV2),  (8) 

and  (Th  (m)  is  the  r.m.s.  height  deviation  determined  from  wind  speed,  p  (m/s),  by 

an  =  O.OOSlp^  (9) 
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In  equation  (8),  /o(.)  is  the  modified  Bessel  function  of  the  first  kind  of  order  zero.  A 
simpler  form  in  terms  of  elementary  functions  has  been  recommended  in  CCIR  Report 
1008-1  [6]: 

Pnfp;  CTh)  ^  ,  (10) 

p.2x  -  2  +  ^(3.2x)2  -  7x  +  9 

where  x  =  1  shows  the  comparison  between  the  exact  and  the  approximate 

po  as  a  function  of  x-  Excellent  agreement  seen  between  the  two  permits  the  use  of  the 
much  simpler  approximate  form  in  actual  numerical  computations.  The  formula  for  po 
makes  physical  sense  only  for  p  >  0.  As  in  the  smooth  earth  case,  we  will  extend  the 
domain  of  definition  of  po  to  include  p  <  0.  We  define  po  for  negative  p  as 


Poi-P]  <^h)  = 


Po{p]  <^h)' 


p  >  0 


We  will  present  the  development  below  for  the  important  practical  case  of  horizontal 
polarization.  Vertical  polarization  can  be  handled  in  a  similar  fashion.  In  the  subsequent 
discussion,  we  will  drop  the  term  CTh  from  the  argument  of  po(p;  (Th)  and  its  presence  will 
be  implied.  For  ranges  far  exceeding  the  heights,  the  plane  waves  arrive  at  very  shallow 
grazing  angles  and  it  is  adequate  to  take  rs(p)  Si  -1  for  horizontal  polarization.  With 
this  approximation,  we  use  (6)  to  write  the  transform  pair  as 

u{x,p)  =  J Pq{p)  f  u{x,z)e'^^  dz - .  f  u{x,z)e  dz  +  w{x,p)  (12a) 

0  y  Poip)  0 

=  =  (-) 

The  additional  function  w{x,p)  must  be  included  to  make  the  pair  defined  in  (12)  consis¬ 
tent.  Of  course,  it  is  not  needed  for  smooth  earth.  Alternately,  one  may  add  a  function 
w{x,  z)  to  (12b)  and  use  (12a)  without  w.  However,  we  prefer  not  to  do  this,  because,  in 
this  case,  the  solution  cannot  be  expressed  in  the  convenient  form  of  (4).  The  quantity 
w{x,p)  is  determined  by  the  requirement  that  —  u.  It  can  be  obtained  as 

.  OO 

w{x,p)  =  -^  j  u[x,q)K{p,q)dq,  (13) 

0 
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where  the  kernel  K{p,q)  =  K{q,p)  is 


\/po{p)po{q) 


Po{p)  -  Po{q)  ,  1  -  Po{p)po{q) 


p-q 


+ 


p  +  q 


For  p  =  q-,  the  kernel  is 


(14) 


K{p,p) 


1 

Po{p) 


Poip)  + 


1  -  pUp) 

2p 


Note  that  (12a)  and  (13)  constitute  a  Fredholm  integral  equation  of  the  second  kind  for 
the  unknown  u  in  terms  of  the  known  function  u.  It  can  be  solved  by  any  one  of  the  several 
techniques  available.  A  simple  technique  to  use  is  the  method  of  successive  approximation 
or  the  so  called  Neumann  series  approach.  In  this  work,  we  only  concentrate  on  some 
simple  techniques  which  are  not  necessarily  efficient  computationally.  Note  that  the  field, 
u,  does  not  satisfy  the  boundary  condition  u(x,0)  =  0  as  it  will  for  the  smooth  earth 
case.  However,  this  does  not  pose  any  conceptual  or  numerical  problems.  For  smooth 
earth,  po  =  1  and,  consequently,  K  =  0.  Also,  observe  that  both  u{x,p)  and  w{x,p)  are 
odd  functions  of  p.  We  may  rewrite  w{x,p)  in  (13)  in  an  operator  form  as 


w{x,p)  =  —iW{u). 


Denoting  the  identity  operator  by  I  and  substituting  in  (12a),  we  arrive  at  the  desired 
forward  transform  as 


u{x,p)  =  =  (I  +  iW)  ^ 


\/po{p)  j  u{x,z)e'^"dz 
0 


CO 

J  u(x,  z)e~'^^  dz 
0 


Defining  the  function  u+(x,z)  by 

u+[x,z)-<^  0  2<0 

and  denoting  its  ordinary  complex  Fourier  transform  as  u+, 

OO 

u+(x,p)  =  J^o(u)  =  j  u+{x,z)e'^'‘  dz, 
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we  write  the  complete  Fourier  transform  pair  for  rough  surface  and  horizontal  polarization 
as 


u(x,p)  =  JF-(w)  =  (2"  + zW)  ^  V^^w+(3:,p)  -  vVo(-p)w+(x, -p) 

u(x,  z)  = 

.yPo{P). 

where  is  the  ordinary  complex  Fourier  inverse 

1 


(15a) 

(156) 


It  remains  to  define  the  inverse  operator  V  =  {I  +  iW)“F  In  computing  the  ordi¬ 
nary  Fourier  transform  pair,  one  would  normally  use  an  A^-point  FFT.  Let  us  assume 
that  the  various  quantities  are  bandlimited  over  — Pmax  ^  P  ^  Pmaxj  S'Od  that  the  trans¬ 
form  is  evaluated  at  p  =  0,  Ap,  2Ap, . . . ,  (iV  -  l)Ap.  Positive  wavenumbers  occur  at 
p  =  Ap,  2Ap, .  •  • ,  (f-  —  l)Ap,  while  negative  wavenumbers  occur  at  (y  -h  l)Ap,  (y  -t- 
2)Ap, . . . ,  (iV  —  1) Ap.  The  value  yAp  corresponds  to  both  -fpmax  and  — Pmax-  Let  us  as¬ 
sume  that  u{x,  ±Pmax)  =  0.  This  can  always  be  achieved  in  practice  by  a  suitable  window 
in  the  transform  domain.  Since  u  is  an  odd  function  of  p,  we  also  have  {t(x,0)  =  0.  The 
same  is  true  of  w.  Consider  the  vectors  u  =  [fi(x,  Ap),  it(x,2Ap), . . .  ,£i(x,  iV/2  -  lAp)]‘ 
and  w  =  [w{x,  Ap),w{x,  2Ap), . . . ,  w{x,  N/2  -  1  Ap)]‘.  One  may  use  trapezoidal  rule  to 
discretize  the  integral  in  (13)  as 


w{x,  mAp)  « 


—iAp 

2x 


N/2-l 

u{x,nAp)K{mAp,nAp). 

n=l 


This  can  be  written  in  a  matrix  form  as 


w  =  —iW  ■  u. 


where  W  of  order  {N/2  -lx  N/2  -  1)  is  the  discrete  version  of  the  continuous  operator 
W  with  elements 

Wmn  =  ^K{mAp,nAp),  m,n  =  1, . .  .,N/2  -  1. 
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The  discrete  form  of  the  continuous  operator,  V,  is  then  P  =  (/  +  W)~\  which  is  a 
matrix  of  order  N/2  -lx  N/2  -  1.  Note  that  W  is  a  full  matrix  and  I  is  an  identity 
matrix  of  size  {N/2  —  lx  N/2  -  1).  In  actual  computations,  we  use  (15a)  to  compute  u 
for  positive  p  and  extend  the  function  over  all  p  using  odd  symmetry.  The  operator  P 
then  acts  on  the  positive  wavenumber  part  of  the  spectrum  occuring  within  the  square 
brackets  of  (15a).  We  give  below  approximations  to  P  of  increasing  orders  of  complexity 
and  accuracy: 


P  =  (/  +  iW)~^  (exact,  0{N^))  (16a) 

/  (zeroth  order,  no  additional  operations)  (165) 

«  /  —  iW  (first  order,  0{N))  (16c) 

«  /  —  iW  —  (second  order,  0{N^))  (16d) 


«  /  -  i0.6438055W  -  0.5936575W^  (least-square,  second  order,  0{N^))  (16e) 

Note  that  the  P  operator  need  be  computed  only  once  using  any  one  of  the  above  ap¬ 
proximations  for  a  given  value  of  wind  speed.  In  subsequent  range  stepping  (equation 
(15a)),  except  for  the  zeroth  order  approximation,  inclusion  of  w  involves  one  additional 
multiplication  of  the  matrix  P  with  the  vector  u  and  is  of  order  0{N'^).  Since  an  FFT 
is  of  order  A'’ln(A^),  the  final  algorithm  time  for  range  stepping  will  be  of  order  O(iV^). 
Diagonalyzing  P  has  the  potential  of  reducing  the  matrix  multiplication  time  to  0{N) 
and  the  overall  algorithm  time  to  0{N\xi{N)).  This  will  be  explored  in  the  future. 

3.  Numerical  Results 

We  validate  our  PE  formulation  first  by  running  a  few  sample  cases  and  comparing  with 
results  from  literature.  A  Hanning  window  with  sequence 

u  \ -I  1  0<n<2|l 

“  \  sin'  ^  <  n  <  f 

is  used  both  in  the  spatial  and  wavenumber  domains.  A  mirror  image  of  h{n)  about 
n  =  0  is  used  for  negative  wavenumbers.  Note  that  the  Hanning  window  forces  a  gradual 


rolloff  to  zero  over  the  last  quarter  of  the  positive  wavenumber  spectrum.  We  label  the 
present  method  as  PEERS,  which  stands  for  Parabolic  Equation  Exact  Rough  Surface. 


The  first  example  we  try  is  that  of  propagation  in  standard  atmosphere  with  M  = 
340  +  0.1 18z,  where  the  height  z  is  in  meters.  The  transmitter  is  at  a  height  of  Zj  =  30 
m,  the  frequency  of  operation  =  3  GHz,  the  maximum  angle  ^max  that  determines  Pmax  is 
1.43°  (=  25  mrad),  maximum  height  including  rolloff  is  z^ax  =  512  m.  The  size  of  EFT 
is  determined  from  z^ax  and  ^max  by 


N>2 


22jnax  sin  0  max 


Ao 


where  Aq  is  the  free-space  wavelength.  The  height  and  wavenumber  increments  are  cal¬ 
culated  from  Nyquist  criterion  as 


Az  = 


2Zn 


N 


and  Ap  = 


27r 

NAz' 


We  choose  the  vertical  increment  Az  =  2  m  and  iV  =  512.  The  range  step  was  chosen 
to  be  200  m.  Figure  2  shows  propagation  factor  (excess  signal  over  the  free-space  case) 
versus  receiver  height  for  at  the  horizontal  range  of  40  km.  Propagation  factor,  PF  in 
dB,  is  defined  as 

PF  =  101og(|u|^xAo), 


where  x  is  the  horizontal  distance  from  transmitter.  A  positive  (negative)  value  of  prop¬ 
agation  factor  implies  gain  (loss)  with  respect  to  propagation  in  free-space.  Although  a 
direct  comparison  is  not  shown  here,  the  lobing  pattern  compares  very  well  with  that 
given  in  [1]. 

The  next  example  we  consider  is  that  of  propagation  in  a  tri-linear  duct  specified  by 


Miz)  = 


340-b0.118z  0<z<135 

499.03  -  1.06z  135  <  z  <  150 

322.33  +  0.118Z  z  >  150. 


where  z  is  in  meters.  The  frequency  of  operation  is  3  GHz,  transmitter  height  is  30  m, 
^max  =  1-43°,  and  Zn,ax  =  512  m.  We  choose  N  =  512,  Az  =  2  m,  and  Ax  =  200  m. 
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Figure  3  shows  propagation  versus  receiver  height  at  a  range  of  40  km.  Once  again  a 
favorable  agreement  with  [1]  was  obtained. 

The  final  example  we  consider  for  smooth  earth  is  propagation  over  a  strong  evap¬ 
oration  duct.  The  frequency  is  10  GHz  and  the  transmitter  height  is  at  25  m.  Re- 
fractivity  profile  for  the  evaporation  duct  is  given  in  Table  1.  We  use  Zmax  1^0  ni, 
^’max  =  1-43°,  ^-2  =  0-293  m,  N  =  1024,  and  Ax  =  200  m.  Figure  4  shows  a  comparison 
of  the  numerical  results  with  those  obtained  by  the  waveguide  model  (MLAYER)  [7]  for  a 
receiver  height  of  25  m.  An  excellent  agreement  is  seen  between  the  two,  thus  validating 
our  PE  algorithm. 

We  now  consider  a  wind-roughened  sea  surface  for  the  evaporation  duct  shown  in 
Table  1.  The  wind  speed  is  10  m/s,  which  corresponds  to  an  r.m.s.  sea  height  of  0.51  m 
according  to  (9).  All  other  parameters  remain  the  same  as  before.  Figure  5  comparses 
the  propagation  factor  versus  range  for  smooth  and  rough  seas.  Clearly,  roughness  has 
resulted  in  the  reduction  of  the  specular  component  as  evidenced  by  the  decreased  ex¬ 
cursions  of  the  signal  in  the  interference  region.  Furthermore,  the  loss  increases  at  larger 
ranges.  Figure  6  shows  the  comparison  of  the  present  numerical  results  with  those  pre¬ 
dicted  by  MLAYER.  Figure  7  shows  comparison  with  TEMPER  developed  by  the  authors 
of  [1]  to  also  include  surface  roughness.  It  is  based  on  first  running  the  PE  on  the  smooth 
earth  case  and  then  estimating  the  grazing  angles  using  one  of  the  spectral  estimation 
procedures.  This  grazing  angles  is  reused  back  into  the  PE  to  modify  the  results  for  the 
rough-earth  case.  The  comparison  is  fair  with  all  the  results  being  within  1.5-2  dB  of 
one  another. 

The  results  shown  in  Figure  6  were  computed  using  the  exact  operator  P  given  in 
(16a).  We  will  now  show  the  effects  of  using  other  approximations  for  P.  In  each 
case  we  will  compare  the  numerical  results  generated  by  the  approximate  operator  with 
those  generated  by  the  exact  operator.  Figure  8  shows  the  results  with  the  zeroth  order 
approximation,  which  is  the  crudest.  Note  the  shift  in  far-out  minimum  relative  to  the 
exact  solution  and  the  overall  disagreement  between  the  two  results.  The  accuracy  can 
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be  dramatically  improved  by  using  the  first  order  approximation.  In  this  case  no  matrix 
inversion  is  required.  Figure  9  shows  the  results  with  the  first  order  approximation. 
There  is  some  disagreement  for  very  long  ranges.  Overall,  solution  by  the  first  order 
operator  compares  favorable  with  the  solution  generated  by  the  exact  operator.  Figure 

10  shows  the  results  with  the  second  order  operator.  It  is  seen  that  the  results  are  almost 
coincident  with  the  exact  results.  Although  the  computation  of  the  exact  and  the  second 
order  operators  are  both  of  the  same  complexity  (of  order  0{N^)),  the  second  order 
operator  can  be  used  when  inversion  of  a  full  matrix  is  not  desired. 

The  last  example  we  consider  is  that  of  a  45.7  m  surface  duct  with  refractivity  profile 
given  by 

J  350  -  0.335Z  0  <  z  <  45.7 

~  \  329.36  +  0.1164Z  z  >  45.7, 

where  z  is  one  again  in  meters.  This  example  was  considered  in  [4].  The  transmitter  is  at 
height  of  25  m,  the  frequency  of  operation  =  10  GHz,  and  the  antenna  is  omnidirectional. 
Propagation  factor  is  calculated  at  a  range  of  200  km  for  various  receiver  heights.  Figure 

11  shows  the  results  for  a  wind  speeds  of  10  m/s  and  20  m/s.  In  both  cases  the  results 
are  compared  with  MLAYER  results.  It  is  seen  that  the  agreement  with  the  latter  is 
very  good,  thereby  validating  our  approach. 

Future  Work 

Our  future  work  will  try  to  remove  any  computational  inefficiencies  in  our  algorithm.  One 
might  formulate  the  problem  in  terms  of  one-sided  sine  and  cosine  transforms  instead  of 
the  full  complex  transforms  defined  in  this  work.  More  efficient  schemes  for  computing 
the  inverse  operator  P  are  certainly  worth  pursuing.  Diagonalizing  the  operator  F  is 
another  are  of  further  investigation.  Extention  to  vertical  polarization  is  straightforward 
and  will  be  looked  at  in  the  future. 
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TABLE  1 


.  REFRACTIVITY  DATA  FOR  EVAPORATION  DUCT 


Height  (m) 

Refractivity  (M  units) 

0.000 

340.00 

0.135 

323.00 

0.223 

321.76 

0.368 

320.53 

0.607 

319.31 

1.000 

318.11 

1.649 

316.94 

2.718 

315.83 

4.482 

7.389 

313.91 

12.182 

313.26 

20.000 

312.99 

313.38 

314.81 

317.99 

148.413 

165.000 

300.000 


324.04 

325.76 

339.745 
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Argument  X 


Figure  1.  Exact  and  CCIR  recommended  formula  for  rough  surface  reduction  factor 
Pq.  The  argument  is  x  = 
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Figiire  3.  Propagation  f8u:tor  versus  receiver  height  at  a  range  of  40  km  in  a  tri- 
linear  duct.  Smooth  earth,  z,  =  30  m,  /  =  3  GHz,  horizontal  polajization  and 
omnidirectional  antenna. 
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factor  versus  range  at  receiver  height  =  25  m  in  a  20  m 
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Figure  7.  Comparison  of  propagation  faetor  versus  range  at  receiver  height  =  25 
m  in  a  20  m  evaporation  duct.  PEERS  versus  TEMPER,  z,  =  25  m,  /i  =  10  m/s, 
/  =  10  GHz,  horizontal  polarization  and  omnidirectional  antenna. 
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Figure  9.  Comparison  of  propagation  factor  versus  range  at  receiver  height  =  25  m 
in  a  20  m  evaporation  duct  by  exact  and  first  order  operators,  z,  =  25  m,  /i  =  10 
m/s,  /  =  10  GHz,  horizontal  polarization  and  omnidirectional  antenna. 
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Figiire  10.  Comparison  of  propagation  factor  versus  range  at  receiver  hei^t  -  5 
m  in  a  20  m  evaporation  duct  by  exact  and  secind  order  operators.  Za  =  25  m, 
u  =  10  m/s,  /  =  10  GHz,  horizontal  polarization  and  omnidirectional  antenna. 
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